%--------------------------------------------------------------------------
%
% computes the leading order (in delta) solution to the basic state and
% then compares to the numerical one
%
%--------------------------------------------------------------------------

function base_composite
close all

p = params;
s = base;

max_v = abs(s.y(end-1,:));
h = s.y(end,:);
t = s.x;

ta = logspace(log10(t(2)), log10(t(end)), 30);
h_comp = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * ta) / (p.beta - 1))))) * (p.beta - 1);

semilogx(t, h, 'k', ta, h_comp, 'r+', 'linewidth',1);
xlim([0, max(t) / 2]);
xlabel('$t$', 'interpreter','latex','fontsize',12);
ylabel('$h(t)$', 'interpreter','latex','fontsize',12);
l = legend('numerical','asymptotic', 'location','southwest');
set(l, 'interpreter','latex','fontsize',11);

v0 = (1 - p.beta) * (1 - 1 ./ h_comp);

psi = 1/6 * p.beta * (1 - p.beta);
psi_0 = psi;
for n = 1:10
    psi = psi + 2 * p.beta * (1 - p.beta) * 1 / n^2 / pi^2 * exp(-n^2 * pi^2 * ta ./ h_comp.^2);
    psi_0 = psi_0 + 2 * p.beta * (1 - p.beta) * (-1)^n / n^2 / pi^2 * exp(-n^2 * pi^2 * ta ./ h_comp.^2);
end


v_comp = v0 + p.delta * (-1/2 * (1 - p.beta - v0) .* (p.beta + v0) .* h_comp + psi);

figure;
loglog(t, max_v, 'k', ta, abs(v_comp), 'r+', 'linewidth',1);
xlabel('$t$', 'interpreter','latex','fontsize',12);
ylabel('$|v(h(t),t)|$', 'interpreter','latex','fontsize',12);
l = legend('numerical','asymptotic', 'location', 'best');
set(l, 'interpreter','latex','fontsize',11);

v_bot = v0 + p.delta * psi_0;

v_mean = (v_comp - v_bot) ./ h_comp;


figure;
semilogx(ta, v_mean);